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A novel algorithm based on the optimized decimation of tensor networks with super-orthogonalization 
(ODTNS) that can be applied to simulate efficiently and accurately not only the thermodynamic but also the 
ground state properties of two-dimensional (2D) quantum lattice models is proposed. By transforming the 2D 
quantum model into a three-dimensional (3D) closed tensor network (TN) comprised of the tensor product 
density operator and a 3D brick-wall TN, the free energy of the system can be calculated with the imaginary 
time evolution, in which the network Tucker decomposition is suggested for the first time to obtain the optimal 
lower-dimensional approximation on the bond space by transforming the TN into a super-orthogonal form. The 
efficiency and accuracy of this algorithm are testified, which are fairly comparable with the quantum Monte 
Carlo calculations. Besides, the present ODTNS scheme can also be applicable to the 2D frustrated quantum 
spin models with nice efficiency. 
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I. INTRODUCTION 

Efficient and accurate numerical methods are very crucial 
to tackle the strongly correlated quantum lattice systems. To 
a large class of intriguing correlated electron and spin mod- 
els, analytical techniques are intractable owing to their ex- 
treme complexity and meanwhile, numerical approaches are 
still challenged by the huge Hilbert space that increases ex- 
ponentially with the lattice size. Two decades ago, numeri- 
cal renormalization group algorithms based on density matrix 
for the ground states 1 and thermodynamic properties 2 of one- 
dimensional (ID) systems were proposed, where the thought- 
ful selection rules were suggested for optimally approximat- 
ing the Hilbert space with an effective subspace. Very re- 
cently, efficient representations with tensor networks as well 
as the corresponding algorithms for the two-dimensional (2D) 
quantum models, for instance, the projected entangle pair state 
(PEPS) 3 , the tree tensor network 4 , the multiscale entangle- 
ment renormalization ansatz state 5 , the infinite PEPS 67 , the 
tensor renormalization group (TRG) 8 ~ 10 , and so on, have been 
suggested. Some of them already gained interesting applica- 
tions (e.g. Refs. GjUl). These algorithms are well testi- 
fied for calculating the ground state properties, while the al- 
gorithms for the thermodynamics of the infinite 2D quantum 
models still need to be developed. 

In this paper, we propose the optimized decimation of ten- 
sor networks with super-orthogonalization (ODTNS) to sim- 
ulate efficiently not only the thermodynamic but also the 
ground state properties of 2D quantum spin lattice models. 
Inspired by the projection method of the ground states of 
2D systems 9 and the linearized TRG method for thermody- 
namic properties of ID systems 13 , we represent the finite tem- 
perature density operator of the 2D quantum model with a 
three-dimensional (3D) closed tensor network (TN) that con- 
sists of the initial tensor product density operator (TPDO) 
and the 3D brick-wall TN for the evolution along the imag- 
inary time direction. The finite temperature properties can be 
obtained by linearly contracting the brick-wall TN with the 



corresponding imaginary time length to get the TPDOll. To 
bound the dimension of the TPDO, we develop the Tucker 
decomposition 14 to the TN and propose the network Tucker 
decomposition (NTD) that transforms a TN into the super- 
orthogonal form so that an optimal lower-dimensional approx- 
imation for the bond space can be reached based on the net- 
work singular value spectrum. We testify the efficiency of the 
ODTNS scheme by calculating the thermodynamic proper- 
ties of the unfrustrated spin- 1/2 Heisenberg antiferromagnet 
on single-layer and bilayer honeycomb lattices, and the ob- 
tained results show the great agreement with quantum Monte 
Carlo (QMC) calculations. We also calculate the thermody- 
namic and magnetic properties of the frustrated bilayer model 
to study the effect of frustration. In what follows, we shall 
present the procedure of the ODTNS algorithm with a 2D 
quantum spin system on a honeycomb lattice as a prototype. 



II. TENSOR NETWORK REPRESENTATION OF THE 
FINITE TEMPERATURE DENSITY OPERATOR 

In accordance with the general definition of the TN's, we 
define a TN as a network consisting of the product of tensors 
(T) and vectors (A), as shown in Fig. Q] Graphically, a point 
with some connected bonds represents a tensor; each bond 
represents an index; a bond that connects two points is called 
geometrical bond which means a shared index by two tensors 
that should be contracted; a bond that only connects one point 
is called the physical bond. We restrict here that one bond 
must connect one or two points. If an index is shared by more 
than two tensors (say n), the restriction would always be ful- 
filled by introducing an nth-order super-diagonal tensor. The 
points on the geometrical bonds represent vectors. A TN with 
no physical bonds is called a closed TN [Fig. G](a)], e.g. a 
TN that denotes the partition function of a classical model; a 
TN with N physical bonds is called an open TN, which con- 
tains d N degrees of freedom (where d is the dimension of one 
physical bond, Fig. [T](b)), e.g. a TN that represents a tensor 
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FIG. 1 : (Color online) (a) A closed square TN in which each tensor 
T has four geometrical bonds that connect each other and no phys- 
ical bond. On each geometrical bond there defines a vector A. (b) 
An open honeycomb TN consisting of two inequivalent tensors T L 
and T R , each of which has three geometrical bonds and one physical 
bond. Three inequivalent vectors A , A 1 ' and A 1 " are defined on three 
inequivalent bonds of the TN, respectively. 




FIG. 2: (Color online) (a) The local evolution operator [/!/., is de- 
composed via an SVD into two gates, G^- g and Gf,j, , each of which 
has two physical bonds (i, i' and j,f, black) and one geometrical 
bond (g, blue); (b) Contract the shared physical bonds among G L and 
G R to get tensors T L and T R ; (c) A TPDO with inverse temperature 
t. Note that the singular value vectors A' " " 1 on each geometrical 
bond are not indicated in (b) and (c) for conciseness. 



product state or a tensor product operator. 

The finite temperature density operator of a 2D system can 
be transformed into an open TN. Suppose that the Hamilto- 
nian can be written as H = Y*ijHij> wnere «y is a local 
Hamiltonian of pairs of spins. The partition function Z is the 
trace of the density matrix p = exp( -J3H) with /3 — l/T the 
inverse temperature and k/j = \. By means of the Trotter- 
Suzuki decomposition 15 , the density operator can be written 
as p st [exp(-T ZijH u )] K+l , where ft = (K + 1)t, and r is 
the infinitesimal imaginary time slice. Define a local evo- 
lution operator t/y = exp(-rHij). Then, the density opera- 
tor can be represented as p =* [{\ uj U U ] K+1 = nf=/ Uij 
where q is the Trotter index. By making a singular value de- 
composition (SVD) on Uyj, = (ij\Uij\i'f) where \ij) stands 
for the direct product basis of spins at site i and /', we have 
t/'/v = Yja G% o A G R ., where A is the singular value vector, 
and G L and G R are two local evolution tensors, each of which 
has two physical bonds (i, i' and j,f, respectively) and one 
geometrical bond (g). For a honeycomb lattice, this step is 
depicted in Figs. |2](a) and (b). 

Next, by contracting the shared bonds among G L and G R 



[Fig. |2](b)], we get 

1 H,gig2g3 ~ Za H-gi 
jk 

T U,gig2g3 = X G 1k,g 2 G kl,gj ' (!) 

Jk 

where g\, gi and g3 are three inequivalent bonds on a honey- 
comb lattice [Fig. |2](c)]. The density operator p at an inverse 
temperature r has the form of a TN as 

p... ff j f ... = Tr G (- ■ -// 2 4f Tj-r^AJjf ' ' -) ' (2) 

in which Ttq is the trace over all contracted geometrical 
bonds, and A 1 , A 11 , A 111 are three inequivalent singular value 
vectors with the initial value A . This gives a TPDO, which is 
an extension of the matrix product density operator 16 and the 
tensor product states. In fact, the TPDO is an open TN com- 
prised of the infinite product of two inequivalent tensors T L 
and T R for two sublattices of the honeycomb lattice as well as 
A 1 , A 11 and A 1 " for three inequivalent bonds [Fig. 0(c)], The 
TPDO at finite temperature consists of two parts: the initial 
TPDO and the 3D brick-wall structure formed by the product 
of evolution tensors. We can contract linearly the evolution 
tensors in pairs into the TPDO along the imaginary time di- 
rection. For example for bond g\, we have 

ik,(g\g\)gigi Z—i ihg'i jk,gig2g3' 
j 

T R = V C R T R (T) 

ik,(gig\)g2g3 l_i ij,g\ jKgxgigi' yJ > 
j 

and meanwhile, we get A ! , = A , A[ . The contractions for 

6 gig\ g\ 81 

bonds g2 and g^ are similar. After certain times of contraction, 
we obtain the TPDO at the corresponding inverse temperature. 
Then by tracing all bonds, we can get the partition function Z 
at finite temperature. During the contraction, as the dimension 
of the geometrical bonds is unavoidably enlarged, an optimal 
approximation is needed to bound the bond dimension. In 
the existing algorithms for truncating the bond, a matrix SVD 
on the matricization of the tensor is used and the states with 
D c (dimension cut-off) largest singular values are preserved. 
We extend the Tucker decomposition to the TN's and suggest 
the NTD to transform the TN into a super-orthogonal form, 
with which the optimal approximation can be obtained with 
the robust network singular value spectrum (NSS). 

III. SUPER-ORTHOGONAL FORM OF TENSOR 
NETWORKS AND NETWORK TUCKER DECOMPOSITION 

In the areas of data compression, image processing, etc., 
Tucker decomposition has been accepted as a convincing 
higher-order generalization of matrix singular value decom- 
position, and its approximation scheme for a single tensor has 
wide and successful applications 17 . It can be written as the 
product of the form 

T iliT< = ^ S JU2 ..,U^U% 2 .-.U^. (4) 

hh-in 
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The tensor S is called the core tensor and lf:\ is the unitary 
matrix. This decomposition is considered as a higher-order 
generalization of the matrix SVD when the core tensor S sat- 
isfies the following two conditions: 

(a) All-orthogonal: 2i 1 / 2 -i a _ 1 i a+1 -4 $ hiT--i*--i* s hir<-< = 
if i a + i' a for any or; 

(b) Ordering: || S ia=l \\>\\ S ia=2 [|> • • ■ >|| S ia=ln ||, where /„ 
is the dimension of the index i n , and the norm of the sub-tensor 

II Si a =k IN Hz, / 2 "i„_i ;„,■•■;„ S hh-ia-\ki„+\S i l i 2 -i c ,-,ki„ +l - 

All-orthogonality requires that each slice (that means fixing 
one index when setting others as one composite index free) of 
the core tensor S is mutually orthogonal with respect to the 
scalar product of matrices. The ordering condition guarantees 
that the norm of each sub-tensors of S does not increase as the 
corresponding index increases, which is similar to the order of 
matrix singular values. Actually, || 5,- a || is the singular values 
of the matrix Mj- la = r,,,^.. .,„...,„ where the composite index 
j = (hh ' ' ' iff-iiff+i • • • in)- 

In the Tucker decomposition, the information of the weight 
is stored in the core tensor, and more specifically, it is the 
norm of each sub-tensor of S . The optimal lower-dimensional 
approximation of a single tensor can thus be obtained by keep- 
ing the space corresponding to the sub-tensors with larger 
norms. Several algorithms of the Tucker decomposition have 
been proposed, such as the higher-order orthogonal iteration 
in which the interplay among all bonds of the tensor is con- 
sidered for the optimal approximation. 

In the following, we extend the definition and the approx- 
imation scheme of the Tucker decomposition for a single 
tensor to a TN. First we define the network reduced matrix 
(NRM) M of bond g, for a (real) tensor T as 



Mgjg. — ^ Tp, gi g2--g i --g„Tp ygigl --g' r - grl (A gi A g2 



(5) 



A-A+i AJ A g, A s'r 



where p = {p\,p%, ■ ■ ■ ,p,„} denotes the composite bond of all 
physical indices because one can always rearrange all physi- 
cal indices into a composite index, g, denotes a geometrical 
bond, and T PiglgT .. g .... gli represents an element of the tensor T. 
The super-orthogonal form of an open TN is defined by two 
conditions: 

(a) Ordering: all A's on geometrical bonds are positive- 
defined, normalized and the elements of each A are in descend- 
ing order. We coin the vectors A's as the network singular 
value spectrum, that is a generalization of the matrix singular 
value spectrum. 

(b) Orthogonality: for any tensor T in the TN and any geo- 
metrical index g, of T, the NRM M is diagonal and equals to 
the square of the corresponding A, say M glg ', = A g 5g ig ' . 

The super-orthogonal conditions, which are non- 
local, require that the matrix A{ pglgr .. g . _ lgM - gn }, gi = 



- />„?lS2'--S,/ L Sl /L g: 



A„. /L. 



A.--A-.. 



A gn (that is analog to the 
singular vectors of matrix SVD) is column orthogonal for any 
;. These conditions are global constraints for the TN, as every 
tensor should satisfy them simultaneously. From the NSS 
which contains the information of the weight distribution 
instead of the core tensor in the Tucker decomposition, the 
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FIG. 3: (Color online) The identical transformation within the NTD 
for the bond gi of an open honeycomb TN. 



optimal low-dimensional approximation of the bond space 
can be obtained. For ID systems, the super-orthogonal 
conditions require that the matrix product states (MPS)'s 
satisfy simultaneously the left and right canonical conditions 
defined in Ref. lfl8ll . which leads to the canonical form of 
MPS. 

Within the suggested NTD, the super-orthogonal form is 
gained by iteratively transforming the TN with identical trans- 
formations on each geometrical bond until the pre-established 
convergence to the super-orthogonal form is reached. For in- 
stance, for bond g\ (Fig. the transformation matrices X L 
and X R for T L and T R (the transformation matrices Y L and Y R 
for bond g2, Z L and Z R for bond g3 are similar) are defined by 



c 

x R ah = 2 Qca(x R r c m u R hL 



(6) 



where t/^ R) contains eigenvectors and Xc (R) contains the 

eigenvalues of the matrix = (A[A{)- 1 M L ^ } . P (Q) con- 

tains the left (right) singular vectors of the intermediate matrix 
W defined by 



w ac = Yj^WAKtft 12 = £ p«$q*. (7) 

b g 

By inserting the identity / = X L {X L y l = X R (X R )~ l into the 
left (right) side of A 1 and transforming T L , T R and A 1 , we have 

Tp, g bc = Tp abc (X L ) a l, ; T R ghc = £ T R ahc (X R ) a l g ; (8) 

a a 

A g 6 gg , = £ X L g XK' = Z P °sW af Q fg ' ■ (9) 

a af 

With new f L , f R and 1' as well as A" and A 1 ", the NRM 
of the bond gi equals to (A 1 ) 2 as the matrices A L p R bc = 

Za'T^tfAfU^fx^^a 112 and P(Q) are both column 
orthogonal and normalized. 
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FIG. 4: (Color online) The convergence of fi, a and £ of open (a) 
honeycomb and (b) square TN's with the increase of the iteration 
step t. [D p \ D g] , /J,,, ■ ■ • , D gn ] is the size of the tensors that form the 
TN, where D gi (i = 1,2, ••■,«) is the dimension of the geometrical 
bonds and D p is the dimension of the physical bond. Each factor is 
obtained from the average of the results of 100 randomly initialized 
TN's. The error of the eigenvalue decomposition <p is about 10 -15 . 

After doing similar transformations on the other two bonds, 
A bears the form of 

-L(R) _ Y ~ m yrpIIrMR), v L(fi) r l/2 
A p,abc ~ /_i 1 p,a'bc A b A c U a'a { X >a 
a' 

= Z C^WWW do) 

a'b'c' 

It can be seen from Eqs. (0 and ( TTOb that the super-orthogonal 
conditions are satisfied when the eigenvalues ^ L(S) in Eq. (O 
are uniformly distributed (i.e., all eigenvalues are equal to 1) 
for all three bonds. The deviation of x UR ^ from uniform dis- 
tribution can be measured by £ = {\x L - *V| + \x R - r V\)/(2£), 
in which | • | means the norm of a vector, *V is a vector with 
all its elements equal to 1 and X is the length of the vector 
X- In addition, we define a factor that measures the conver- 
gence of A's at the fth iteration by //(f) = TiS=i n,m(\^ (' _ 
3) - A s (t)\/\A s (t)\)/3, and a factor cr = (cr 1 + a*)/!, where 
cr L{R) = Y,ab \M L a b R) - U S a ) 2 Sab\ measures TN's deviation from 
the super-orthogonal form according directly to the super- 
orthogonal conditions. 

To testify the robustness of the super-diagonal form and 
the efficiency of NTD, we randomly initialize the inequiva- 
lent tensors and vectors (according to Gaussian distribution 
N{Q, 1)) that form the infinite open honeycomb and square 
TN, and calculate the factors fi, cr and f with different iter- 
ation steps. The value of each factor is the average of the re- 
sults of 100 randomly initialized TN's. Fig. |4]shows that the 
NTD can transform a randomly initialized TN into a super- 
orthogonal form very efficiently. It is found that ( and cr 
decay exponentially to about 10~ 14 ~ 10~ 15 (where the er- 
ror of the eigenvalue decomposition itself q is about 10~ 15 ) 
within 300 steps. Meanwhile, yu converges to 10~ 14 ~ 10~ 15 , 
which justifies the good convergence of the three inequivalent 
/l's. We may see that for a certain TN, the three factors share 
one same super-orthogonalization ratio £ and in general, £ be- 
comes larger when we increase the space of the physical bond 
and fix the space of the geometrical bonds. 

The computational cost to transform a TN into the super- 



orthogonal form with NTD is mainly from the eigenvalue (sin- 
gular value) decompositions of D g x D g matrices, which is 
about 0(tD 3 g ) with t the transformation steps and D g the di- 
mension of the geometrical bond. In the ODTNS scheme, the 
TPDO converges to the super-orthogonal form with cr < 10~ 8 
only with f ~ 10 steps, because for a small r, the evolution is 
nearly identical. 



IV. THE FREE ENERGY 

After each time that the evolution tensors are contracted 
into the TPDO, we super-orthogonalize the TPDO and obtain 
the optimal approximation of the enlarged geometrical bonds 
as well as the A s 's (S = I, II, HI). By collecting the nor- 
malization factor = -^Tjg with q the Trotter step, the 

free energy per site can be obtained with and 7 19 that is the 
contraction of the TPDO by 




The thermodynamic quantities of the 2D quantum lattice sys- 
tems can be obtained from the free energy f(J3). 

What is more, the ground state properties can also be ob- 
tained with the ODTNS scheme. When one takes K — » oo and 
r — » 0, the ground state energy per site eo has a simple form 
of 

e = lim lim — In ]~[ r s . (12) 



V. THERMODYNAMICS OF SPIN-1/2 
ANTIFERROMAGNET ON A HONEYCOMB LATTICE 

To judge the efficiency and accuracy of the ODTNS algo- 
rithm, let us consider the spin- 1/2 antiferromagnet on an in- 
finite honeycomb lattice with H u = 6(SfS*. + + 
where 5 measures the anisotropy of spin interactions. In fol- 
lowing calculations, // is kept smaller than 10~ 820 , and the lat- 
tice size for QMC calculations is 64 x 64. 

Fig. [5] shows the energy difference AE — \E - Eqmc\ be- 
tween the results obtained by ODTNS and QMC calculations 
at different inverse temperature fi in the absence of the mag- 
netic field, where r = 10~ 2 and D c = 30. We find that, when 
6 = 0.5, there exists a thermodynamic phase transition, and 
the energy difference is about 10~ 4 ~ 10~ 5 at both high and 
low p. Near the critical point the difference is relatively high 
but is still smaller than 0.003. When 5-1, the system is gap- 
less, and the energy difference is about 10~ 3 at both high and 
low ft, and near the crossover point, the difference is also rel- 
atively high but still remains around 10~ 2 . These results show 
that the precision of ODTNS scheme is comparable with that 
of QMC. 

We investigated AE versus the dimension cut-off D c near 
the critical (crossover) ft, as shown in Fig. |6l where r = 10" 2 . 
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FIG. 5: (Color online) The energy difference AE = \E — Eq M c \ be- 
tween the ODTNS and QMC calculations and the QMC error at dif- 
ferent inverse temperature /? for the spin- 1/2 Heisenberg antiferro- 
magnet on a honeycomb lattice, (a) shows the result at S = 0.5 when 
the system suffers a thermodynamic phase transition; (b) shows the 
result at 5 = 1 when the system is gapless and a thermodynamic 
phase transition is forbidden by Mermin-Wagner theorem. We set 
t = 10~ 2 ,o- < 10~ 8 and D c = 30. 
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FIG. 6: (Color online) The energy difference AE = \E — Eq M c \ be- 
tween ODTNS and QMC with the dimension cut-off D c near the crit- 
ical (crossover) point for (a) 6 = 0.5 and (b) 5 = 1. It can be seen that 
AE becomes smaller as D c is increased. The QMC error is around 

io- 5 . 



It is observed that the energy difference becomes smaller 
when D c is increased. When /3 is away from the critical 
(crossover) point, we uncovered that different D c gives errors 
within 10 4 . We also checked AE for different t, and dis- 
closed the (Trotter) errors are within 10~ 4 . 

The specific heat as a function of /3 is calculated by C = 
-/3 2 dE/d/3, as shown in Fig. |7]for 6 = 0.5. A divergent peak 
at a critical temperature T c is observed, which indicates that 
a phase transition occurs between a paramagnetic phase and 
an antiferromagnetic phase at T c . Such a phase transition is 
also confirmed with the result of the staggered magnetization, 
shown in the inset of Fig. Q The QMC results are also in- 
cluded for a comparison. One may see that both results from 
the ODTNS and QMC calculations agree quite well, showing 
again the efficiency and accuracy of the present method. In ad- 
dition, the present ODTNS algorithm can be directly applied 
to the 2D frustrated quantum spin models. 
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FIG. 7: (Color online) The inverse temperature /3 dependence of the 
specific heat at S = 0.5 and h s = for the spin- 1/2 anisotropic 
Heisenberg antiferromagnet on a honeycomb lattice. The QMC re- 
sult with the error around 10~ 3 is included for a comparison. The in- 
set shows the staggered magnetization at the magnetic field h s = 0.01 
and 0.05, and the QMC error is around 10~ 5 . 



VI. THERMODYNAMICS OF A SPIN-1/2 FRUSTRATED 
BILAYER HONEYCOMB HEISENBERG MODEL 



We apply the ODTNS algorithm to explore the spin- 1/2 
Heisenberg model on a bilayer honeycomb lattice with the 



Hamiltonian H 

X«) , Mb) 



2</7> Hlj, 



where H u = + Sff 



(W" + Hf>)/3, = JitlSwiSfS 3 } + S]S]) + S'Sy the 

anisotropic Heisenberg antiferromagnet on each layer and 
fff' b) = Jad^MS-Sj + + the interlayer cou- 

pling [see the inset of Fig. [8] (a) for the layout of the model], 
5i t 2 and 6 a j, measure the corresponding anisotropy of nearest 
neighbor spin interactions. We take J\ — J 2 — 1 as energy 
scale. 

When both J a and Jt, are positive, the couplings are all anti- 
ferromagnetic, and the system has no frustration. The energies 
obtained by the ODTNS algorithm at J' — J a — Jt, — 1 and 
3 are shown in Fig. H](a), which are in good agreement with 
QMC results (where the QMC error is within 10~ 5 ). 

When J a and J/, take different signs, the system becomes 
frustrated and the QMC simulations fail because of suffering 
from the negative sign problem. Fig. |8](b) shows the results 
of energy and specific heat at J' — J a — -Jf, = 1 and 3. 
When /' = J a = 1, a second-order phase transition is found 
at f3 c = 2.75(5); when J' = 3, no phase transition is observed 
from the specific heat. Notice that the frustration reaches the 
maximum at /' = 3 in the Ising limit for this present model. 
Fig. [9] presents the sublattice magnetization per site m s at 
J' — 1 and 3 where the frustration exists. At J' = 1, the 
couplings within each layer are dominant. The system is in 
the antiferromagnetic phase at low temperature and there ex- 
ists a thermal phase transition from the antiferromagnetic to 
paramagnetic phase. At J' = 3 when the strong frustration 
is present, m s is around 10~ 3 , indicating the absence of mag- 
netic long range order at all temperature. These calculations 
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FIG. 8: (Color online) The inverse temperature {S dependence of the 
energy for the spin- 1/2 Heisenberg model on a bilayer honeycomb 
lattice for (a) J' = J a = Jt = 1 and 3, and (b) ]' = ]„ = -Ji, = 1 and 
3. The inset of (a) shows the bilayer structure of the system and the 
inset of (b) shows the specific heat C of the 2D frustrated quantum 
spin system. 
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FIG. 9: (Color online) The sublattice magnetization per site m s at 
J' = J a = 1 and 3 (inset) for the spin- 1/2 frustrated Heisenberg 
model on a bilayer honeycomb lattice, where J a IJi, = -1. 



show that the ODTNS algorithm is capable of studying the 2D 
frustrated quantum spin systems. 



VII. SUMMARY 

In summary, a novel algorithm based on the ODTNS 
scheme for the 2D quantum spin lattice models is proposed. 
By mapping the 2D quantum model into a 3D TN, we sug- 
gest the NTD to obtain the optimal approximation of the bond 
space by transforming the TPDO into the super-orthogonal 
form, that leads to an efficient and accurate calculation of 
the free energy as well as other observables in the 2D quan- 
tum systems. We testify the efficiency and accuracy of the 
present algorithm by studying the thermodynamics of a spin- 
1/2 Heisenberg antiferromagnet on a honeycomb lattice, and 
compare the results with those of the QMC. It is shown that 
the precision of the ODTNS algorithm is comparable with 
that of the QMC as both results agree very well. In addition, 
we find that the present algorithm can also be applied to ex- 
plore the 2D frustrated quantum spin models without suffering 
from a negative-sign problem. It is expected that the present 
ODTNS scheme could also be extended to 2D correlated elec- 
tron systems. 
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Note that without the iteration in the NTD, a is on the same order 
of t for the imaginary time evolution, which could lead to the 
results unstable near the critical (crossover) point. If one reduces 
r to decrease cr without making iterations in the NTD, the cost 



would be very high because the extremely large evolution steps 
are needed to reach a certain temperature. 



